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I.  INTRODUCTION 

A.  BACKGROUND 

Modeling  has  become  a  widely  accepted  way  to  analyze  complex  engineering 
problems.  One  of  the  most  powerful  tools  for  this  is  MATLAB™,  a  computer  programing 
language  specifically  geared  towards  numerical  modeling  of  engineering  problems.  It  has 
a  wide  variety  of  uses  and  can  help  explain  engineering  principles  which  would  be 
difficult  or  too  expensive  to  actually  conduct  in  a  laboratory. 

B.  PURPOSE 

This  thesis  is  an  extension  of  work  done  by  LCDR  Thomas  A.  Green  of  providing 
computer  algorithms  for  topics  in  the  text,  Fundamentals  of  Acoustics,  Third  Edition, 
John  Wiley  &  Sons,  Inc.,  by  Kinsler,  Frey,  Coppens,  and  Sanders.  The  programs  are  to 
provide  an  interactive  tool  for  the  student  to  use  in  broadening  their  understanding  of 
material  in  the  text. 

C.  REQUIREMENTS 

Although  these  programs  can  be  used  without  reference  to  text  material,  they 
will  be  most  useful  when  used  in  conjunction  with  the  sections  of  KFCS  which  apply. 
The  student  will  be  able  to  more  quickly  grasp  the  relationships  between  the  input  and 
output  of  the  program  when  the  text  is  used  in  concert  with  the  computer  as  the  programs 
are  designed  to  further  explain  the  text  material. 

The  programs  in  this  thesis  will  run  on  any  machine  which  supports  MATLAB™, 
Professional  Version  4.2  or  later.  Although  the  student  version  of  MATLAB™  can  be 
used  in  most  cases,  some  of  the  algorithms  use  large  vectors  and  matrices  which  cannot 
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be  handled  by  the  student  version  of  MATLAB™. 

To  obtain  more  information  about  MATLAB™,  contact  The  Math  Works,  Inc., 
University  Sales  Department,  24  Prime  Part  Way,  Natick,  Massachusetts  01760-1500,  or 
call  them  at  (508)  653-1415,  or  go  to  their  World  Wide  Web  site  at 
http://www.mathworks.com. 

A  copy  of  the  programs  in  this  thesis  can  be  obtained  by  sending  a  formatted  disk 

to: 

Professor  James  V.  Sanders 
Physics  Department  (Code  PH) 
Naval  Postgraduate  School 
833  Dyer  Road,  Room  203 
Monterey,  CA  93943-51 17 


EL  DAMAGED  ARRAYS 

The  program  BKNARRAY.M  displays  plots  of  beam  pattern  for  arrays  of  any 
size.  The  most  useful  feature  of  this  program  is  it's  ability  to  illustrate  how  the 
performance  of  an  array  declines  as  elements  within  it  fail.  This  program  supports 
instruction  in  concepts  discussed  in  Chapter  8  of  KFCS. 
A.         CONCEPTS 

The  first  assumption  is  that  the  distance  to  the  receiver  is  much  greater  than  the 
length  of  the  array.  It  is  generally  accepted  that  a  difference  of  one  order  of  magnitude 
satisfies  this  requirement .  Additionally,  the  program  assumes  that  the  elements  radiate 
spherical  waves  with  the  same  amplitude  and  phase.  Although  this  program  was 
developed  as  an  array  of  sources,  by  the  principle  of  reciprocity,  the  beam  patterns  can 
also  be  viewed  as  patterns  of  reception. 

Section  8. 13  of  KFCS  discusses  the  simple  line  array.  Because  of  our  assumption 
that  the  elements  emit  waves  of  the  same  amplitude  and  phase,  the  form  of  the  pressure 
field  of  each  source  is 


A     Kot-kr) 

~*  (2.1) 


Where  r;  is  the  distance  from  the  ith  element  to  the  position  in  the  far  field.  Because  the 


beam  pattern  is  normalized  to  the  pressure  on  the  axis  of  the  main  lobe,  the  amplitude  term 
can  be  ignored  and  the  term  which  contributes  to  the  beam  pattern  of  an  unsteered  array  is 

e  ~]kr'  (2.2) 

which  is  basically  the  phase  contribution  of  each  element.  The  total  pressure  in  the  far  field 
is  then  simply  the  sum  of  the  pressure  contribution  of  all  elements  in  the  array  or 

£  « '""'  (2.3) 

r'=l 

Arrays  are  normally  designed  to  have  a  single  major  lobe  as  narrow  as  possible. 
The  criterion  to  achieve  this  is 

d=^JT  (24) 

where  d  is  the  array  spacing.  This  program  uses  the  optimal  array  spacing  to  obtain  a 
single  major  lobe  so  frequency  is  not  a  user  input. 

Finally,  the  plots  are  referenced  to  the  acoustic  pressure  emitted  by  a  single 
element  and  the  same  dB  scale  is  maintained  for  each  set  of  plots. 
B.         PROGRAM  DESCRIPTION 

The  program  has  two  user-controllable  parameters.  The  first  is  the  number  of 
elements  in  the  array;  this  must  be  a  positive  number.  If  the  user  is  only  interested  in  the 


beam  pattern  for  a  specific  number  of  elements  then  this  is  the  only  input.  To  see  how  the 
array  behaves  as  elements  fail  they  can  input  the  total  number  of  elements  they  would  like 
to  see  fail.  The  program  plots  both  the  beam  pattern  for  a  fully  functioning  array  and  the 
beam  pattern  as  each  additional  element  fails.  The  failed  elements  are  chosen  at  random  by 
the  computer.  Functioning  elements  in  the  array  are  denoted  by  a  green  "o"  and  failed 
elements  are  denoted  by  a  red  "x".  These  are  displayed  above  the  array's  beam  pattern. 
Figure  2. 1  is  a  sample  of  the  plots  generated  by  BKNARRAY.M.  On  the  lefthand 


ARRAY  GAIN 
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0      180 


270  270 

Beampattern  for  1 0  Elements  Beampattern  for  1 0  Elements 

Element  3  has  Failed 


Figure  2.1  Beam  Pattern  for  a  10  Element  Array 


side  of  the  figure  is  a  plot  of  the  beam  pattern  of  a  ten  element  array  and  on  the  right  is 

the  plot  of  the  same  array  with  the  third  element  not  functioning.  BKNARRAY.M  uses  a 

modified  version  of  the  MATLAB™  program  POLAR.M  which  is  included  in  the 

appendix. 

C.        ALGORITHM  bknarray.m 

function  bknarray(elements,failures) 

%  BKNARRAY  is  a  function  used  to  display  the  beampattern  of  an  n-element  array.  The 

%  number  of  elements  is  chosen  by  the  user.  In  addition  the  user  can  disable  any  number 

%  of  elements  in  the  array  to  see  how  this  affects  the  beampattern.  If  it  is  desired  to  see 

%  just  the  fully  functioning  array  pattern  then  the  user  must  input  "0"  as  the  number  of 

%  failed  elements.  The  user  is  not  able  to  select  which  elements  they  want  to  disable, 

%  only  the  total  number.  The  function  randomly  selects  which  elements  to  disable  and 

%  displays  the  array  format  below  the  array  pattern.  The  dB  reference  is  the  gain  above 

%  what  a  single  element  would  see. 

global  radius  %  radius  is  declared  as  a  global  variable  in  order  to  allow  the 

%  function  POLAR3  to  plot  all  beampatterns  of  an  array  on  the 

%same  scale, 
global  r3db 
% 

%  This  is  the  user  protection  portion  of  the  program. 


if  nargin=  1 

failures=0; 

if  (elements  <  1) 

errorCNumber  of  elements  must  be  more  than  0.') 
end 
elseif  nargin  =  2 

if  (elements  <  1) 

errorCNumber  of  elements  must  be  more  than  zero.1) 
end 
if(failures<0) 

errorCNumber  of  failures  must  be  zero  or  greater.') 
end 
if  (failures  >  elements) 

errorCBKNARRAY  requires  failures  <  or  =  elements.') 
end 
else 

errorCBKNARRAY  requires  at  least  one  input  argument.') 
end 


o 


frequency=250;  %  in  hertz 


soundspd=  1 500;  %  in  meters/sec 

k=(2* pi* frequency )/soundspd;  %  wave  number 

lambda=soundspd/frequency;  %  wavelength 

arrayspace=(lambda*(elements- 1  ))/(elements);         %  spacing  of  array  elements 

range=(arrayspace*(elements- 1  )*500);  %  distance  in  meters 

theta=0:(2*pi/180):(2*pi);  %  determines  number  of  analysis  points 

failed_element=[]; 

element=ones(  1  ,elements); 

plot_position=l; 

radius=  1 0*log  1 0(elements); 

if  elements— 1 

r3db=radius-3; 
else  r3db=0; 
end 
element_position=([0:  (elements- 1 )] .  *arrayspace)-. . . 

((elements/2  )*  array  space)+(arrayspace/2 ); 
element_position=element_position' ; 
[THETAJPOS]=meshgrid(theta,element_position); 
%  elementrange  computes  the  range  to  each  analysis  position  from 
%  each  array  element 


element_range=(((range*sin(THETA)).  A2)+. . . 

(((range*cos(THETA))+POS)A2)).A(.5); 
pressure=(exp(-j.*element_range.*k));  %  computes  pressure  contribution 

%  of  each  element 
if  failures— 0 

ele_mat=  1  elements; 
for  ind=l:failures; 

choice=ceil(rand*length(ele_mat)); 
failed_element=[failed_elementele_mat(choice)]; 
ele_mat(choice)=[] ; 
end 
end 

for  count=0:  failures, 
ifcount~=0 

pressure(failed_element(count),:)=zeros(size(pressure(  1 ,:))); 
element(failed_element(count))=0; 
end 
if(elements=l) 

beampattem=  1 0*log  1 0(abs(pressure)); 
else 


beampattern=  1  O*logl  0(abs(sum(pressure))); 

one_gain:=((beampattern>0)|(beampattern=0)); 

beampattern=beampattern.  *one_gain; 

%  sums  pressure  contribution  of  each  array  element 

%  to  obtain  beampattern 
end 

figure(ceil((count+ 1  )/2)) 
if  plot_position=3 

plotjposition=l; 
end 

subplot(  1 ,2,plot_position) 
polar3(theta,beampattern) 
title('ARRAY  GAIN1) 
y=1.6*radius; 
hold  on 
for  ind=l:elements 

x=2*(ind-(elements/2)-.5); 

if  element(ind)=l 
plot(x,y,'go') 

elseifelement(ind)=0 
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plot(x,y,'rx*) 

end 
end 

hold  off 

y=  1 . 1  *radius*sin(45  *pi/ 1 80); 
x=y; 

text(x,y,[num2str(round(radius)),'dB']) 
ifelements>2 

y=l.l*r3db*sin(30*pi/180); 

x=l.l*r3db*cos(30*pi/180); 

textCx^nim^sntround^db)),'  dB']) 
end 
text((-radius),(- 1 . 5  *radius),. . . 

['Beampattern  for  ^nurr^strCelements),'  Elements'],. 

'verticalalignment','bottom') 
ifcount~=0 
text((-radius),(- 1  7*radius),. . . 

['Element  ',num2str(failed_element(count))v . . 

'  has  Failed'],'verticalalignment','bottom') 
end 
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plot_position=plot _position+ 1 ; 
end 
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m.  RESONANCE  IN  PIPES 

The  program  FPIPE.M  computes  the  acoustic  resistance  and  reactance  of  a 
flanged,  open  ended  pipe.  Plots  are  generated  which  can  be  used  to  determine  resonant 
and  antiresonant  frequencies  for  the  pipe.  This  program  supports  instruction  in  concepts 
developed  in  Chapter  9,  Section  2  of  KFCS. 
A.         CONCEPTS 

To  simplify  the  analysis  of  resonance  within  the  flanged  pipe  the  analysis  is 
restricted  to  frequencies  such  that  ka«l,  where  k  is  the  wave  number  and  a  is  the  radius 
of  the  pipe  so  only  plane  waves  will  propagate  in  the  pipe.  The  frequencies  of  resonance 
for  any  mechanical  system  occur  at  the  points  where  the  imaginary  component  of  input 
mechanical  impedance,  the  input  mechanical  reactance,  becomes  zero.  If  the  input 
mechanical  resistance  is  non-zero,  these  zero  points  can  also  define  antiresonance. 
Antiresonance  and  resonance  can  be  separated  by  calculating  the  real  component  of  input 
resistance  at  these  zero  points.  Resonant  points  are  those  with  small  values  of  resistance 
and  antiresonant  points  have  large  values  of  resistance.  The  equation  for  input  mechanical 
impedance  of  a  pipe  with  characteristic  impedance  p0c  and  cross  section  S  is 

—OL+jtankL 


'-       P4 (3-D 


pncS  Z  , 

pQcS 
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where  Z^  for  a  flanged  pipe  is  given  by 


ZmL       l/Lo,8 


(ka)  y-—ka  (3  2) 

p0c5    2  3tc  v      7 

By  substituting  (3.2)  into  (3.1)  and  solving  for  Zm0/(p0cS)  the  resonant  and  antiresonant 
frequencies  of  the  pipe  can  be  found. 
B.         PROGRAM  DESCRIPTION 

The  program  has  three  arguments  for  user  input.  These  are  pipe  length,  pipe 
diameter,  and  sound  speed.  All  three  of  these  arguments  are  required  to  be  entered  for  the 
program  to  function.  Pipe  length  should  be  entered  in  meters  while  pipe  diameter  should 
be  in  centimeters  and  sound  speed  should  be  in  meters/second. 

The  program  determines  a  maximum  frequency  by  using  the  requirement  that 
wavelength  be  much  larger  than  the  radius  of  the  pipe.  This  requirement  generally  means 
that  the  wavelength  be  an  order  of  magnitude  larger  than  the  pipe  radius  and  FPIPE.M 
uses  this  criterion  to  set  the  frequency  limit. 

The  values  for  reactance  and  impedance  have  been  normalized  by  the 
characteristic  impedance  of  the  fluid.  This  does  not  affect  the  frequencies  of  resonance 
or  antiresonance  so  they  are  accurate  for  any  fluid  medium  with  the  stated  sound  speed. 

Figure  3.1  is  a  plot  of  the  input  reactance  for  a  1.5  meter  long  pipe  which  is  6 
centimeters  in  diameter  in  water  with  a  sound  speed  of  1500  m/s.  The  plot  has  been 
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magnified  using  the  "zoom"  command  to  show  the  zero  crossings.  Figure  3.2  shows  a 
plot  of  the  input  resistance  of  the  same  pipe.  The  points  where  the  reactance  and  the 
resistance  become  zero  are  the  resonant  points.  Antiresonance  occurs  where  the 
resistance  becomes  a  maximum. 


Flanged  Pipe  1 .5  Meters  Long,  Sound  Speed  1500  m/s 
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Figure  3.1  Input  Reactance  for  a  Flanged  Pipe. 
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Flanged  Pipe  1 .5  Meters  Long.  Sound  Speed  1500  m/s 
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Figure  3.2  Input  Resistance  for  a  Flanged  Pipe. 

C.         ALGORITHM  fpipcm 

function  fpipe(plength,pdiam,sndspd) 

%  FPIPE(pipe  length,pipe  diameter,sound  speed)  This  program  determines  the  load 
%  resistance  and  reactance  for  a  flanged  pipe.  Because  it  is  assumed  that  the  wavelength 
%  is  much  larger  than  the  transverse  dimension  of  the  pipe,  the  program  only  plots 
%  frequencies  whose  wavelength  is  an  order  of  magnitude  greater  than  the  pipe  width 
%  chosen.  The  resistance  and  reactance  information  is  plotted  against  frequency  in  two 
%  separate  plots.  This  information  can  then  be  used  to  find  values  of  resonance  and 
%  anti-resonance  for  the  pipe,  "plength"  should  be  in  meters,  "pdiam"  should  be  in 
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%  centimeters,  and  "sndspd"  should  be  meters/sec. 
%This  is  the  user  protection  portion  of  the  program, 
if  nargin~=3 

error(TPIPE  requires  three  input  arguments.') 
end 

chk_input=min([plength,pdiam,sndspd]); 
if(chk_input<=0) 

error('All  arguments  of  FPIPE  must  be  positive.'); 
end 

%  This  is  the  beginning  of  the  computations. 
pipe_radius=pdiam/(2  *  1 00); 

max_freq=round(sndspd/(pipe_radius*10));  %  long  wavelength  requirement 
freqs=  1 : .  5 :  maxrreq; 
k=(2 .  *pi.  *freqs)./sndspd; 

%  This  section  finds  the  input  mechanical  impedance  assuming  that  the  pipe  is  flanged. 
Zml=((((k.*pipe_radius)A2)./2)+(j.*(8/(3*pi)).*(k.*pipe_radius))); 
ZmO=(Zml+(j.*tan(k.*plength)))./(l+(j*(Zml.*(tan(k.*plength))))); 
Rm0=real(Zm0); 
Xm0=imag(Zm0); 
zero_ref=zeros(size(freqs)); 
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figure(l) 

plot(freqs,(RmO)); 

ylabel('Input  Resistance') 

xlabelCFrequency  (Hz)') 

title([Tlanged  Pipe  ',num2str(plength),'  Meters  Long,  Sound  Speed ',. 

num2str(sndspd),'  m/s']) 
figure(2) 

plot(freqs,zero_ref) 
hold  on 

plot(freqs,(XmO)); 
hold  off 

ylabel('Input  Reactance') 
xlabel(Trequency  (Hz)') 
title(  [Flanged  Pipe  ',num2str(plength),'  Meters  Long,  Sound  Speed ',. 

num2str(sndspd),'  m/s']) 
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IV.  POWER  RADIATION 

The  program  PEPEPWR.M  computes  the  power  radiated  from  open  ended  pipes  of 
different  sizes  for  both  flanged  and  unflanged  terminations.  This  program  supports 
concepts  developed  in  Chapter  9,  Section  3  of  KFCS. 
A.         CONCEPTS 

For  a  pipe  of  radius  a  terminated  in  a  flange  the  transmission  coefficient  is 


T_=- 


2{kaf 


n 


\+k*af?+(±?(kaf  (41) 

2  371 


Similarly  the  transmission  coefficient  for  an  unflanged  pipe  is 


T=- 


(tar)2 


+±(ka)2]2+(0.6ka)2  (42) 


These  equations  are  valid  when  ka«  1  so  that  the  denominator  is  nearly  one  and  the  two 
equations  become 

T^2(kaf  (4.3) 

for  the  flanged  pipe  and 

T^(ka)2  (4.4) 
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for  the  unflanged  pipe.  PEPEPWR.M  computes  all  of  these  transmission  coefficients  and 
plots  them  to  allow  the  user  to  see  the  difference  that  this  assumption  makes. 
B.         PROGRAM  DESCRIPTION 

The  program  uses  only  a  single  argument,  pipe  diameter,  which  should  be  in 
centimeters.  Figure  4. 1  is  a  plot  of  the  transmission  coefficient  using  both  equations  for  a 
flanged  pipe.  Figure  4.2  is  the  same  type  of  plot  for  an  unflanged  pipe. 


Coefficients  for 
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Figure  4.1  6cm  Diameter  Flanged  Pipe  Transmission  Coefficient 


Notice  that  in  both  the  flanged  and  unflanged  case  the  estimated  transmission 
coefficient  increases  more  rapidly  than  the  exact  coefficient.  This  is  because  we  have 
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ignored  the  contribution  of  ka  in  the  denominator  and  as  frequency  increases  this  term 
becomes  more  important. 


Coefficients  for  an  Unflanged  Pipe  6  Cm  in  Diameter 
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Figure  4.2  6cm  Diameter  Unflanged  Pipe  Transmission  Coefficient 

Figure  4.3  shows  the  power  transmission  coefficient  for  a  flanged  pipe  1 
centimeter  in  diameter.  Notice  that  the  shape  of  the  plot  is  very  similar  to  that  of  Figure 
4. 1  which  is  a  pipe  6  centimeters  in  diameter.  It  is  evident  that  the  transmission 
coefficient  of  the  1  cm  pipe  is  much  lower  than  that  of  the  6  cm  pipe  at  the  same 
frequency.  This  illustrates  the  idea  that  small  acoustic  devices  are  poor  transmitters  of 
relatively  long  wavelength  signals. 
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Coefficients  for  a  Flanged  Pipe  1  CM  in  Diameter 
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Figure  4.3  1cm  Diameter  Flanged  Pipe  Transmission  Coefficient 


Finally  by  comparing  Figures  4. 1  and  4.2  the  flanged  pipe's  estimated  coefficient 
is  double  that  of  the  unflanged  pipe's  coefficient  as  is  expected.  Also  as  expected  the 
actual  transmission  coefficient  of  the  flanged  pipe  is  much  less  than  twice  the  actual 
coefficient  of  the  unflanged  pipe. 
C.        ALGORITHM  pipepwr.m 
function  pipepwr(diam) 

%  PIPEPWR  computes  the  power  transmission  coefficient  for  small  diameter  pipes 
%  which  are  flanged  and  unflanged.  The  upper  curve  of  each  plot  is  the  exact  solution 
%  when  the  assumption  is  only  that  the  wavelength  is  much  larger  than  the  diameter  of 


22 


%  the  pipe.  The  lower  curve  is  the  solution  when  it  is  assumed  that  the  wavenumber 
%  times  the  radius  of  the  pipe  is  much  less  than  one.  The  argument  "diam"  is  the 
%  diameter  of  the  pipe.  Values  for  diam  should  be  in  units  of  centimeters. 
%  This  is  the  user  protection  portion  of  the  program, 
if  nargin~=l 

error(PIPEPWR  requires  only  one  input  argument.1) 
end 

chk_input=min(diam); 
ifchk_input<=0 

error('All  values  for  radius  must  be  positive.') 
end 

%  This  begins  the  body  of  the  program. 
radius=diam/(2*100); 

max_freq=343/(radius*  10);    %  Low  frequency  limit 
freqs=l  :(round(max  freq)); 
k=2*pi.*freqs./343; 
ka2=(k.*radius)A2; 

T_flange=(2.*(ka2))./(((l+(-5*ka2)).A2>+(((8/(3*pi))A2).*(ka2))); 
T_fl_est=2*ka2; 
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T_unflange=(ka2)./(((l+(.25*ka2)).A2)+(-36*ka2)); 

T_unfl_est=ka2; 

figure(l) 

plot(freqs,T_flange,'y-') 

hold  on 

plot(freqs,T_fl_est/r:') 

hold  off 

xlabel(Trequency') 

ylabel(Power  Transmission  Coefficient') 

title(['Coefficients  for  a  Flanged  Pipe ',... 

num2str(diam),'  CM  in  Diameter1]) 
legendCExact  Coeff '/Est  Coeff ') 
figure(2) 

plot(freqs,T_unflange,'y-') 
hold  on 

plot(freqs,T_unfi_est,'r:') 
hold  off 

xlabel(Trequency') 
ylabel(Power  Transmission  Coefficient') 
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title( ['Coefficients  for  an  Unflanged  Pipe  ',... 

num2str(diam),'  Cm  in  Diameter']) 
legend(*Exact  Coeff*,'Est  Coeff* ) 
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V.  STANDING  WAVES 

The  program  STNDWAVE.M  uses  the  single  argument  of  standing  wave  ratio 
(SWR)  to  compute  the  load  impedance  over  a  range  of  nodal  position  values.  The  data  is 
then  plotted  and  the  load  impedance  for  a  given  SWR  and  node  position  can  be  read  from 
the  graph.  This  program  supports  instruction  in  concepts  discussed  in  KFCS  Chapter  9, 
Section  4. 
A.         CONCEPTS 

A  load  impedance  will  cause  a  reflection  of  an  acoustic  wave  from  the  load.  The 
ratio  of  the  reflected  pressure  wave  to  the  incident  pressure  wave  in  a  pipe  of  length  L  and 
cross  sectional  area  S  can  be  determined  from  the  SWR  by  the  following  equation 

B    SWR- I 


A    SWR  +  1 


(5.1) 


By  substituting  (5.1)  into 


-        1  ^e<» 
— =^—  (5  2) 


the  complex  load  impedance  normalized  by  the  characteristic  impedance  p0c  can  be 
determined.  If  0,  the  phase  angle  between  the  reflected  and  incident  waves,  is  known  the 
equation  determining  the  phase  angle  is 
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Q=2k(L-xl)-n 


(5.3) 


where  x,  is  the  position  of  the  node  closest  to  the  load.  The  argument  (L-x,)  can  be  scaled 
so  that  it  is  measured  as  a  fraction  of  a  wavelength.  This  will  allow  cancellation  with  the 
k  term  and  make  the  determination  of  8  independent  of  frequency. 
B.         PROGRAM  DESCRIPTION 

The  program  uses  only  the  single  user  input  of  SWR.  The  results  obtained  from 
this  program  apply  to  any  fluid  medium,  any  sound  speed,  and  any  length  pipe.  It  also 
applies  to  any  frequency  or  pipe  diameter  as  long  as  the  long  wavelength  criterion  (X»a) 
is  met.  Figure  5. 1  shows  the  load  resistance  of  a  pipe  with  a  SWR  of  4.  The  resistive 

Load  Resistance  of  a  Pipe  with  a  SWR  of  4 
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Figure  5.1  Load  Resistance  with  SWR=4 
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load  reaches  a  peak  when  the  node  is  one  quarter  wavelength  from  the  end.  This  peak 
resistive  load  corresponds  to  the  point  where  the  maximum  energy  is  being  transmitted  to 
the  medium.  Figure  5.2  is  a  plot  of  the  load  reactance  with  a  standing  wave  ratio  of  4. 


Load  Reactance  of  a  Pipe  with  a  SWR  of  4 
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Figure  5.2  Load  Reactance  with  S WR=4 


Notice  that  the  reactance  passes  through  zero  at  the  same  quarter  wave  point  where  the 
resistance  reached  its  peak.  Figures  5.3  and  5.4  are  plots  of  load  resistance  and 
impedance  for  a  SWR  of  2.  Notice  that  as  the  SWR  goes  down  the  load  resistance  and 
impedance  also  go  down. 
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Load  Resistance  of  a  Pipe  with  a  SWR  of  2 
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Figure  53  Load  Resistance  for  a  SWR=2 


Load  Reactance  of  a  Pipe  with  a  SWR  of  2 
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Figure  5.4  Load  Reactance  for  a  SWR=2 
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C.        ALGORITHM  stndwavcm 

function  stndwave(SWR) 

%  STNDWAVE  computes  and  plots  the  load  impedance  for  a  pipe  with  a  given 

%  measured  standing  wave  ratio  over  a  range  of  L-x  distances.  L-x  is  the  distance  from 

%  the  load  end  of  the  pipe  to  the  first  acoustic  node.  This  distance  is  measured  as  a 

%  fraction  of  a  wave  length.  Since  there  are  two  nodes  per  wavelength,  it  is  only 

%  necessary  to  plot  up  to  one  half  of  a  wavelength.  SWR  is  a  unitless  ratio  of  the 

%  acoustic  pressure  at  a  node  to  that  at  an  antinode. 

%  This  is  the  user  protection  portion  of  the  program. 

if  nargin~=l 

error('STNDWAVE  requires  only  one  input  argument.') 
end 
ifSWR<l 

error(' Values  for  SWR  must  be  equal  to  or  greater  than  1.') 
end 

%  This  begins  the  body  of  the  program. 
off_dist=0:.005:.5; 
press_ratio=(S  WR- 1  )/(S  WR+ 1 ); 
theta=((4*pi).*off_dist)-pi; 
Z=(  1 +(press_ratio.  *exp(j .  *theta)))./(  1  -(press_ratio.  *exp(j .  *theta))); 
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x=[0,.5]; 

y=[0,0]; 

figure(l) 

plot(off_dist,real(Z)) 

xlabel(T)istance  from  Termination  in  Fraction  of  a  Wavelength') 

ylabelCLoad  Resistance') 

titletfLoad  Resistance  of  a  Pipe  with  a  SWR  of  ',num2str(SWR)]) 

figure(2) 

plot(off_dist,imag(Z)) 

hold  on 

plot(x,y) 

hold  off 

xlabel(Distance  from  Termination  in  Fraction  of  Wavelength') 

ylabelCLoad  Reactance') 

title(['Load  Reactance  of  a  Pipe  with  a  SWR  of  ',num2str(SWR)]) 
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VI.  COMBINED  DRIVER-PIPE  SYSTEMS 

The  program  SYSRES.M  determines  resonance  frequencies  of  combined  driver 
and  pipe  assemblies.  It  plots  system  reactance  against  a  scale  normalized  by  kL  which 
allows  the  user  to  determine  resonance  conditions.  This  program  supports  concepts 
discussed  in  Chapter  9,  Section  6  of  KFCS. 
A.         CONCEPTS 

When  combining  two  or  more  elements  into  a  system  it  is  important  to  understand 
how  the  components  interact.  Earlier  chapters  of  KFCS  addressed  the  resonance  of  drivers 
and  sections  of  Chapter  9  address  the  resonance  of  pipes.  When  these  are  combined  the 
system  will  resonate  when  the  reactance  of  the  driver  is  equal  to  the  reactance  of  the  pipe. 
Equation  (6.1)  shows  the  relationship  of  driver  impedance  to  pipe  impedance. 

(coskL)(sinkL)  ,  f  _  b 

sm2kL+(aLfcos2kL~a       ~kL  (61) 

Where  a  is  the  ratio  of  the  moving  mass  m  of  the  driver  to  the  mass  of  the  fluid  in  the  pipe. 


m 

a=  (6.2) 


SpJL 
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b  is  the  ratio  of  the  stiffness  s  of  the  suspension  of  the  driver  to  the  stiffness  of  the  fluid  in 
the  pipe. 


b=- 


sL 


SPoc2 


(6.3) 


a  is  the  absorption  coefficient  of  the  pipe. 
B.         PROGRAM  DESCRIPTION 

The  program  requires  two  inputs  from  the  user,  the  values  of  a  and  b.  Figure  6. 1  is 
a  plot  of  the  reactances  for  a  light,  flexible  driver  and  Figure  6.2  is  a  plot  for  a  heavy,  stiff 
driver. 


15 


10 


Resonances  of  a  Combined  Pipe  -  Driver  System 
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Figure  6.1  Light  Flexible  Driver. 
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Resonances  of  a  Combined  Pipe  -  Dnver  System 
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Figure  6.2  Heavy  Stiff  Driver 
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The  following  values  were  used  for  Figure  6. 1;  a=04,  and  b=2.67.  For  Figure  6.2  the 

values  were;  a=25,  and  b=32.  Where  two  lines  cross  corresponds  to  a  system  resonance. 

It  must  be  remembered  that  the  values  you  pick  for  k  need  to  satisfy  the  long  wavelength 

criteria  for  the  analysis  to  be  valid. 

C.        ALGORITHM  sysres.m 

function  sysres(a,b) 

%  SYSRES  computes  and  plots  the  resonant  points  of  a  combined  driver  and  pipe 

%  assembly,  "a"  is  the  ratio  of  the  mass  of  the  driver  to  the  mass  of  the  fluid  in  the  tube. 

%  "b"  is  the  ratio  of  the  stiffness  of  the  driver  to  the  stiffness  of  the  fluid  in  the  tube.  The 

%  tube  is  closed  at  one  end. 
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%  This  is  the  user  protection  portion  of  the  program, 
if  nargin~=2 

error('SYSRES  requires  two  input  arguments. ') 
end 

chk_input=Tnin(  [a,b] ); 
if  chk_input<0 

error('Arguments  for  SYSRES  must  zero  or  greater.') 
end 

%  This  begins  the  body  of  the  program. 
domain=5*pi; 
dom_int=domain/200; 
kL=dom_int:dom_int:domain; 
alpha_L=01; 
plot_dom=5; 
plot_dom_int=5/200; 

KL=plot_dom_int:plot_dom_int:plot_dom; 
react_driver=(a.  *kL)-(b./kL); 

react_pipe=(cos(kL).*sin(kL))./(((sin(kL)).A2)+((alpha_L).A2).*. 
(cos(kL))); 
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max_react=l  .2*max(react_pipe); 

if  max_react>15  max_react=15; 

end 

zero_pty=[0,0]; 

zero_ptx=[0,plot_dom] ; 

figure(l) 

plot(KL,react_dri  ver,': ') 

hold  on 

plot(KL,react_pipe,'-') 

plot(zero_ptx,zero_pty,V) 

hold  off 

axis([0,plot_dom,-maxjreact,maxjreact]) 

titleCResonances  of  a  Combined  Pipe  -  Driver  System') 

xlabelCValue  of  kL  Scaled  by  Pi') 

ylabelCLoad  and  Driver  Impedance') 

legend(*Driver  '/Pipe',-1) 
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VH.  RECTANGULAR  CAVITY 

The  program  RECTCAV.M  displays  the  nodal  pattern  of  a  rectangular  cavity.  The 
user  inputs  the  dimensions  of  the  cavity  and  specifies  the  desired  mode.  The  program 
computes  the  pattern  along  with  the  frequency  to  excite  that  mode.  This  program  supports 
instruction  in  concepts  discussed  in  Chapter  9,  Section  7  of  KFCS. 
A.         CONCEPTS 

The  boundary  conditions  of  a  rigid-walled  rectangular  cavity  require  that  the 
pressure  of  a  normal  mode  be  a  combination  of  cosines  so  the  acoustic  pressure  is  a 
maximum  at  the  walls  and  solutions  for  the  wave  equation  will  be  of  the  form 

/'faB^JcosW(coV)(cosy^"'  V-  V 

with  frequencies  given  by 


(xilmn = Jcos2(Itz/Lx)  +cos2(imz/L  ) +cos2(mz/L7) 

B.         PROGRAM  DESCRIPTION 

The  program  RECTCAV.M  requires  six  input  arguments  from  the  user.  These  are 
the  lengths  of  the  three  sides  of  the  cavity  and  the  number  of  nodes  in  each  of  the  three 
directions.  The  program  then  determines  the  resonant  frequency  of  the  given  mode  and 
plots  the  nodal  pattern  in  three  dimensions.  The  user  can  use  the  "view"  command  to 
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change  which  side  of  the  cavity  is  being  viewed.  By  changing  the  axis  limits  the  user  can 

look  inside  the  box  to  see  the  amplitude  at  a  specific  location.  The  program  plots  in  color 

and  as  such  will  not  reproduce  in  black  and  white.  In  the  plot,  red  is  maximum  amplitude 

or  antinodes.  A  red  area  surrounded  by  blue  and  green  is  180  degrees  out  of  phase  with  a 

red  area  surrounded  by  yellow  and  orange.  Light  blue  corresponds  to  nodes. 

C.         ALGORITHM  rectcav.m 

function  rectcav(length,  width,  height,  1,  m,  n) 

%  RECTCAV  computes  the  resonant  frequency  of  a  rectangular  cavity  with  a  given 

%  mode.  It  then  plots  the  nodal  pattern  and  displays  the  resonant  frequency. 

%  This  is  the  user  protection  section. 

ifnargin~=6 

error( 'RECTCAV  requires  six  arguments.') 
end 

chk_input=min([length,width,height]); 
if  (chk_input<0)|(chk_input=0) 

error('Dimensions  of  cavity  must  be  greater  than  0.'); 
end 

chk_input=min([l,m,n]); 
if(chk_input<0) 

error('Values  for  l,m,  and  n  must  be  0  or  greater.'); 
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end 

if  ((l=0)&(m=0)&(n=0)) 

error('At  least  one  value  of  l,m,  or  n  must  be  greater  than  zero.1); 
end 

%  Start  of  computations. 
no_steps=30; 
sndspd=343; 

freq=((((l*pi/length)A2)+((m*pi/width)A2)+((n*pi/height)A2))A5); 
freq=rreq*(sndspd/(2*pi)); 
colormap(hsv); 

max_length=max([length,width,height]); 

x=linspace(0,(l*pi),no_steps);  %  Computes  arguments  of  cosine 

y=linspace(0,(m*pi),no_steps); 
z=linspace(0,(n*pi),no_steps); 
if  (1=0) 

x=zeros(  1  ,no_steps); 
end 
if(m=0) 

y=zeros(  1  ,no_steps); 
end 
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if(n=0) 

z=zeros(  1  ,no_steps); 
end 

Lx=linspace(0,length,no_steps);  Ly=linspace(0,width,no_steps); 
Lz=linspace(0,height,no_steps); 
[X,  Y]=meshgrid(x,y);  xysurf=(cos(X)).  *(cos(Y)); 
co  efr=ones(size(xy  surf)) ; 
figure(l) 
hold  on 
for  ind=l:(size((Lz),2)) 

Z=coeff.*(Lz(ind)); 

xyzsurf=xysurf.  *(cos(z(ind))); 

mesh(Lx,Ly,Z,xyzsurf) 
end 

axis([0,max_length,0,max_length,0,max_length]) 
view(3) 

title(['Rectangular  Cavity  with  Resonant  Frequency  of  ,num2str(round(freq)),'  Hz']) 
xlabel('Length  of  Cavity1) 
ylabel('Width  of  Cavity') 
zlabel('Height  of  Cavity') 


42 


Vm.  RECTANGULAR  WAVEGUIDES 

The  program  GRPSPD.M  computes  the  group  speed  and  phase  speed  in  a 
rectangular,  rigid-walled  waveguide.  This  program  supports  instruction  in  concepts 
developed  in  KFCS  Chapter  9,  Section  8. 
A.         CONCEPTS 

The  modes  found  in  a  waveguide  for  a  specific  drive  frequency  are  given  by 

(v>lcf=kl+kt  (8.1) 

where  0)  is  the  driving  frequency,  k^,  is  the  wave  number  of  the  propagating  wave,  and  kta 
is  the  wave  number  of  the  lm  mode.  If  oo/c  <  k^  then  the  lm  mode  is  evanescent.  The  wave 
number  of  a  propagating  mode  is 


^(u)/c)2<  (8.2) 


k 


The  equation  for  the  group  speed  is  given  by  equation  (8.3). 


c^HcdJg))2  (8.3) 

So  energy  travels  down  the  waveguide  at  a  speed  less  than  the  speed  of  sound. 

The  phase  speed,  which  is  how  fast  a  surface  of  constant  phase  will  travel  in  the 
waveguide,  is 
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Cp  m — rv  (8-4) 


Phase  information  travels  down  the  waveguide  at  a  speed  greater  than  the  speed  of 
sound.  In  fact  the  phase  speed  becomes  infinite  when  the  driving  frequency  equals  the 
cutoff  frequency. 
B.        PROGRAM  DESCRIPTION 

The  program  requires  the  user  to  provide  values  of  waveguide  height  and  width 
and  mode  number.  The  height  and  width  should  be  in  meters  and  at  least  one  of  the  mode 
numbers  should  be  greater  than  zero.  The  program  plots  the  straight  line  at  a  ratio  of  one 
which  corresponds  to  plane  wave  propagation  when  both  phase  speed  and  group  speed 
are  equal  to  the  speed  of  sound.  The  program  also  plots  the  mode  input  by  the  user  and  a 
mode  which  is  two  more  in  height  and  one  more  in  width.  Figure  8. 1  is  a  plot  of  the  (1,1) 
mode  and  the  (3,2)  mode  for  a  6  cm  by  6  cm  waveguide.  As  the  driven  frequency 
approaches  infinity  the  phase  and  group  speeds  for  all  modes  are  equal  to  the  sound 
speed.  The  (0,0)  mode  is  a  plane  wave  for  all  drive  frequencies  so  its  phase  speed  and 
group  speed  are  both  equal  to  the  speed  of  sound  for  all  frequencies.  This  graph  can  be 
used  to  determine  group  speed  and  phase  speed  for  any  mode  in  any  fluid  medium. 
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Group  and  Phase  Speeds  Relative  to  Sound  Speed 
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Figure  8.1  Group  and  Phase  Speeds  of  a  Rectangular  Cavity  6cm  by  6  cm. 


C.         ALGORITHM  grpspd.m 

function  grpspd(height,width,l,m) 

%  GRPSPD  plots  the  relation  of  group  and  phase  speed  to  drive  frequency.  The  user  is 

%  required  to  provide  values  for  waveguide  height,  waveguide  width,  and  mode 

%  numbers.  Height  and  width  should  be  in  meters  and  at  least  one  mode  number  should 

%  be  greater  than  0. 

%  This  is  the  user  protection  portion  of  the  program. 

ifnargin~=4 

error('GRPSPD  requires  four  input  arguments.') 
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end 

if  min(  [height,  width]  )<=0 

errorCHeight  and  width  of  waveguide  must  be  positive.') 
end 
ifmin([l,m])<0 

errorCMode  numbers  must  be  zero  or  greater.1) 
end 
if(l=O)&(m=0) 

error('At  least  one  of  the  mode  numbers  must  be  greater  than  0.') 
end 

%  This  begins  the  body  of  the  program. 
12=1+2; 
m2=m+l; 

f_co=round((343/(2*pi))*((l*pi/height)A2  +  (m*pi/width)A2)A  5); 
f_co2=ceil((343/(2*pi))*((12*pi/height)A2  +  (m2*piAvidth)A2)A5); 
f_end=10*f_co2; 

freqs  1  =logspace(log  1 0(f_co),log  1 0(fend),  1 00); 
freqs2=logspace(logl0(f_co2),logl0(f_end),100); 
freqsO=logspace(  1 , log  1 O(fend),  1 00); 
wl=2*pi.*freqsl; 
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w2=2*pi.*freqs2; 

wcol=2*pi*f_co; 

wco2=2*pi*f_co2; 

groupO=ones(size(freqsO)); 

groupl=real((l-((wcol./wl).A2))A.5); 

group2=real((l-((wco2./w2)A2)).A.5); 

phase  1=1.  /group  1; 

phase2= 1  ./group2; 

figure(l) 

semilogx(freqs  1  ,group  1 ,'-') 

hold  on 

semilogx(freqs  1  ,phase  1 ,':') 

semilogx(freqs2,group2,'-') 

semilogx(f^eqs2,phase2,,:,) 

semilogx(freqsO,groupO,'w-') 

hold  off 

axis([10,f_end,0,2]) 

titleCGroup  and  Phase  Speeds  Relative  to  Sound  Speed') 

xlabel(T)rive  Frequency') 

ylabelCRatio  of  Group  or  Phase  Speed  to  Sound  Speed') 
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textC^co/lO^.S^nunLZstitD/^numlstrlmX'Mode']) 
text(f_co2*2,0.5,[num2str(12),,,,,num2str(m2),,Mode']) 
legend('Group  Speed',Phase  Speed',- 1) 
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EX.  HELMHOLTZ  RESONATOR 

The  Helmholtz  resonator  is  an  example  of  a  lumped  acoustic  element.  Where  all 
dimensions  are  much  smaller  than  a  wavelength.  The  two  programs  HELMRES 1  M  and 
HELMRES2.M  determine  resonant  frequencies  of  Helmholtz  resonators  over  a  range  of 
neck  lengths  and  a  range  of  neck  opening  sizes  respectively.  These  programs  support 
instruction  in  concepts  developed  in  KFCS  Chapter  10,  Section  2. 
A.         CONCEPTS 

A  Helmholtz  resonator  is  a  cavity  with  an  opening.  The  frequency  of  its  resonance 
can  be  determined  by  establishing  the  impedance  characteristics  of  the  cavity.  The 
reactance  of  the  resonator  is  controlled  by  a  mass-type  element  and  a  stiffness-like 
element.  For  the  Helmholtz  resonator  the  mass  is  provided  by  the  fluid  in  the  neck. 
Because  the  neck  is  open  end  effects  must  be  accounted  for  so  the  equation  of  effective 
mass  is  given  by:  m=p0SL'  where  L'  for  a  flanged  neck  as:  L'=L+1.7a  ,  a  being  the 
diameter  of  the  neck.  For  an  unflanged  neck  the  equation  is:  L'=L+1.5a.  The  stiffness  is 
determined  by  the  volume  V  of  the  cavity,  as  the  fluid  in  the  cavity  acts  like  a  spring  as  it  is 
compressed.  The  stiffness  is  given  by  (9. 1) 

*=Poc2y  (9.1) 


where  S  is  the  area  of  the  neck.  So  now  the  reactance  can  be  determined  by  substituting 
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into  the  equation 

Xm=(wm-s/u>)  (92) 

and  solving  this  when  reactance  equals  zero  gives 


G)0=c 


(9.3) 


\L'V 


This  relationship  holds  for  any  cavity  as  long  as  the  long  wavelength  criterion  is  met. 
B.         PROGRAM  DESCRIPTION 

The  programs  HELMRES1.M  and  HELMRES2.M  both  compute  resonant 
frequencies  for  cavities  and  require  three  inputs  from  the  user.  For  HELMRES1  M,  the 
user  is  required  to  input  values  of  sound  speed,  neck  diameter,  and  cavity  volume  while 
HELMRES2.M  requires  sound  speed,  neck  length,  and  cavity  volume.  Both  programs 
plot  resonances  for  flanged  and  unflanged  necks  but  only  the  flanged  neck  results  are 
shown  here.  Figure  9. 1  shows  the  resonant  frequencies  over  a  range  of  neck  lengths  and 
Figure  9.2  shows  the  frequencies  over  a  range  of  neck  diameters. 
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Resonant  Frequencies  of  Helmholtz  Resonators 
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Figure  9.1  Resonances  of  Different  Neck  Lengths 
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Resonant  Frequencies  of  Helmholtz  Resonators 
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Figure  9.2  Resonances  for  Different  Neck  Openings 
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C.         ALGORITHM  helm  res  l.m 

function  helmres  1  (sndspd,hole_diam,volume) 

%  HELMRES  1  computes  resonant  frequencies  of  a  Helmholtz  Resonator  as  a  function  of 

%  neck  length,  "sndspd"  should  be  in  units  of  meters/sec,  "holediam"  should  be  in 

%  centimeters,  and  "volume"  should  be  in  cubic  centimeters. 

%  This  is  the  user  protection  portion  of  the  program. 

if  nargin~=3 

errorCHELMRESl  requires  three  input  arguments.') 
end 

chk_input=min(  [sndspd,hole_diam,  volume] ); 
ifchk_input<=0 

error('All  arguments  of  HELMRES  1  must  be  positive  numbers.') 
end 

%  This  starts  the  body  of  the  program. 
hole_diam=hole_diam/ 1 00 ; 
volume=volume/ 1  e6; 
radius=hole_diam/2; 
L=0:.01:l; 

flangeL=L+(  1 .7*radius); 
unflangeL=L+(  1 . 5  *radius); 


52 


S=pi*(radiusA2); 

flange_arg=((S./(flangeL.*volume)).A.5); 

unflange_arg=((S./(unflangeL.*volume)).A.5); 

flange  freq=flange_arg.*(sndspd/(2*pi)); 

unflange_freq=unflange_arg.*(sndspd/(2*pi)); 

figure(l) 

plot(  flange_freq,L,' W) 

xlabelCFrequency  of  Flanged  Neck  (Hz)') 

ylabelCLength  of  Neck  (meters)') 

titleCResonant  Frequencies  of  Helmholtz  Resonators') 

figure(2) 

plot(unflange _freq,L,'w') 

xlabelCFrequency  of  Unflanged  Neck  (Hz)') 

ylabelCLength  of  Neck  (meters)') 

titleCResonant  Frequencies  of  Helmholtz  Resonators') 

D.        ALGORITHM  helmres2.m 

function  helmres2(sndspd,neck_length,  volume) 

%  HELMRES2  computes  resonant  frequencies  of  a  Helmholtz  Resonator  as  a  function  of 

%  neck  opening  size,  "sndspd"  should  be  in  units  of  meters/sec,  "necklength"  should  be 

%  in  centimeters,  and  "volume"  should  be  in  cubic  centimeters. 
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%  This  is  the  user  protection  portion  of  the  program, 
if  nargin— 3 

error(,HELMRES2  requires  three  input  arguments.') 
end 

chk_input=min(  [sndspd,neck_length,  volume] ); 
ifchk_input<=0 

errorCAll  arguments  of  HELMRES2  must  be  positive  numbers.') 
end 

%  This  starts  the  body  of  the  program. 

volume=volume/le6;  %  converts  volume  to  cubic  meters 

radius=.001:.001:.025; 

Radius=.l:.  1:2.5; 

L=neck_length/100;  %  converts  neck  length  to  meters 

flangeL=L+(  1.7.  *radius);        %  These  compute  the  equivalent 

unflangeL=L+(1.5.*radius);    %  neck  lengths. 

S=pi.*(radiusA2); 

flange_arg=((S./(flangeL.*volume)).A.5); 
unflange_arg=((S./(unflangeL.*volume)).A.5); 
flange  _freq==flange_arg.*(sndspd/(2*pi)); 
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unflange _freq=iinflange_arg.*(sndspd/(2*pi)); 

figure(l) 

plot(flange_freq,Radius,'w') 

xlabelCFrequency  of  Flanged  Neck  (Hz)') 

ylabelCRadius  of  Neck  Opening  (centimeters)') 

titleCResonant  Frequencies  of  Helmholtz  Resonators') 

figure(2) 

plot(unflange  freq^adius/w1) 

xlabel(Trequency  of  Unflanged  Neck  (Hz)') 

ylabelCRadius  of  Neck  Opening  (centimeters)') 

titleCResonant  Frequencies  of  Helmholtz  Resonators') 
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X.  RESONANT  BUBBLES 

A  small  bubble  in  water  can  be  viewed  as  a  lumped  acoustic  element.  The  two 
programs  explore  the  resonant  frequencies  and  the  acoustic  losses  of  small  air  bubbles  in 
water.  The  program  only  analyzes  to  a  depth  of  1 00  meters  as  there  are  not  many  air 
bubbles  of  any  consequence  below  this  depth.  These  programs  support  instruction  in 
concepts  explored  in  Chapter  10,  Section  3  of  KFCS. 
A.         CONCEPTS 

As  in  the  case  of  the  Helmholtz  resonator,  the  bubble  can  be  treated  as  a  lumped 
acoustic  element.  Therefore  its  resonant  frequency  is  controlled  by  the  same  mechanical 
properties  which  lets  us  determine  that  the  resonant  frequency  is 


co0= 


1 


-LJ.  (io.i) 


where  a  is  the  radius  of  the  bubble,  y  is  the  ratio  of  specific  heat  of  gas  in  the  bubble,  T0  is 
the  ambient  pressure,  and  p0  is  the  density  of  the  water. 

These  small  bubbles  can  be  important  in  determining  the  losses  of  a  transmitted 
acoustic  wave  because  they  absorb,  reflect,  and  reradiate  the  sound  energy.  The  total 
effect  of  energy  losses  due  to  the  bubble  is  described  by  its  extinction  cross  section.  It  is 
defined  as  the  ratio  of  energy  absorbed  and  scattered  to  the  intensity  of  the  incident 
energy.  This  can  be  broken  down  into  a  scattering  cross  section  and  an  absorption  cross 
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section.  The  ratio  of  the  scattering  cross  section  to  the  extinction  cross  section  gives  us  an 
idea  of  whether  scattering  or  absorption  is  the  more  dominant  loss  mechanism  and  can  be 
found  by  using  the  equation 


— =fy*2  (10.2) 

o 


where  kois  the  wave  number  of  the  resonant  frequency  in  water  and  the  quality  factor  Q  is 
given  by 

2=  (10.3) 


^a+l^xlO"4^ 


B.         PROGRAM  DESCRIPTION 
1.  bubble.m 

This  program  is  used  to  determine  the  frequency  at  which  a  bubble  of  a  specified 
diameter  will  resonate.  The  diameter  of  the  bubble  should  be  entered  by  the  user  in 
millimeters.  The  depths  analyzed  range  from  0  to  100  meters  depth.  It  is  important  to 
realize  that  the  program  determines  the  resonance  for  a  bubble  with  a  constant  diameter 
whatever  the  depth.  However,  a  bubble  which  starts  at  a  certain  size  at  100  meters  depth 
will  grow  larger  as  it  rises  and  so  this  will  also  contribute  to  a  change  in  the  resonant 
frequency.  Figure  10.1  is  a  plot  of  the  resonant  frequencies  of  two  bubbles,  one  with  a 
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Resonant  Frequency  of  an  Air  Bubble  in  Seawater 
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Resonant  Frequency  (Hz) 

Figure  10.1  Resonant  frequencies  of  an  air  bubble  2  mm  in 
diameter  and  one  of  3mm  in  diameter. 

diameter  of  2  millimeters  and  one  with  a  diameter  of  3  millimeters.  From  equation  (10. 1) 
it  can  be  seen  that  as  the  bubble  diameter  decreases  the  resonant  frequency  will  increase, 
so  the  2mm  bubble  is  the  righthand  curve  and  the  3mm  bubble  is  the  lefthand  curve.  The 
user  can  input  as  many  different  sizes  of  bubbles  as  desired  in  a  single  chart  by  using 
vector  notation  to  define  the  argument  of  bubble  diameter. 


2. 


bubble2.m 


This  program  determines  the  ratio  of  scattering  cross  section  to  extinction  cross 
section  at  resonant  frequency  which  allows  us  to  determine  which  loss  mechanism  is 
dominant.  Figure  10.2  shows  a  plot  of  this  ratio  for  the  same  two  bubbles  of  2  mm  and 
3  mm  in  diameter.  Again  the  2  mm  bubble  is  on  the  left  and  the  plot  of  the  3  mm  bubble  is 
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Scattering  Cross  Section/Extinction  Cross  Section 
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Figure  10.2  Ratio  of  scattering  cross  section  to  extinction  cross 
section  for  a  2mm  diameter  bubble  and  a  3mm  diameter  bubble. 


on  the  right.  It  can  be  seen  from  the  graph  that  as  depth  increases  the  scattering  losses 
become  more  important  than  the  absorption  losses.  Again  the  user  can  plot  as  many 
bubble  sizes  on  the  same  plot  as  desired  by  using  vector  notation. 
C.         ALGORITHM  bubble.m 
function  bubble(radius) 

%  BUBBLE  computes  the  resonant  frequency  of  air  bubbles  within  the  water  column 
%  down  to  a  depth  of  100  meters.  Any  number  of  values  for  radius  may  be  entered  by 
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%  using  vector  notation,  i.e.  BUBBLE([vall,val2,...,valN]).  This  allows  the  user  to 

%  compare  air  bubbles  of  different  size  on  the  same  graph.  Values  for  radius  should  be  in 

%  units  of  millimeters. 

%  This  is  the  user  protection  portion  of  the  program. 

if  nargin~=l 

error(To  enter  more  than  one  value  for  radius  use  vector  notation.') 
end 

chk_input=min(radius); 
if  chk_input<=0 

error('All  values  for  radius  must  be  positive.') 
end 

%  This  begins  the  body  of  the  program. 
depth=l:100; 
radius=radius./1000; 
spec_heat=1.01; 
rho=1026; 

pressures. 013e5+((1026*9.81).*depth); 
[P,R]=meshgrid(pressure,radius); 
omega=((3*spec_heat.*P./rho).A5)./R; 
freq=omega/(2*pi); 
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figure(l) 

plot(freq,depth,V) 

axis('ij') 

xlabel('Resonant  Frequency  (Hz)') 

ylabel('Depth  Below  Surface  (meters)') 

title('Resonant  Frequency  of  an  Air  Bubble  in  Seawater') 

D.         ALGORITHM  bubble2.m 

function  bubble2(radius) 

%  BUBBLE2  computes  the  ratio  of  the  scattering  cross  section  to  the  extinction  cross 

%  section  for  different  size  air  bubbles  in  the  water  column  down  to  a  depth  of  1 00 

%  meters.  Any  number  of  values  for  radius  may  be  entered  by  using  vector  notation. 

%  i.e.  BUBBLE2([vall,val2,...,valN]).  This  allows  the  user  to  compare  air  bubbles  of 

%  different  size  on  the  same  graph.  Values  for  radius  should  be  in  units  of  millimeters. 

%  This  is  the  user  protection  portion  of  the  program. 

if  nargin~=l 

error(To  enter  more  than  one  value  for  radius  use  vector  notation.') 
end 

chk_input=min(radius); 
if  chk_input<=0 

error('All  values  for  radius  must  be  positive.') 
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end 

%  This  begins  the  body  of  the  program. 

depth=l:100; 

radius=radius./1000; 

spec_heat=1.01; 

rho=1026; 

sndspd=1500; 

pressure=1.013e5+((1026*9.81).*depth); 

coefY=ones(size(pressure)); 

[C,R]=meshgrid(coeff,radius); 

Radius=C.*R; 

[P,R]=meshgrid(pressure,radius); 

omega=((3*spec_heat.*P./rho).A5)./R; 

kO=omega/sndspd; 

Q=l./((kO.*Radius)+((1.6e-4).*(omega.A5))); 

ratio=kO.  *Radius.  *Q; 

figure(l) 

plot(ratio,  depth,  V) 

axis('ij') 

xlabel('Ratio  of  Cross  Sections') 


63 


ylabel('Depth  Below  Surface  (meters)') 

title('Scattering  Cross  Section/Extinction  Cross  Section') 
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XL  REFLECTED  WAVES  IN  PIPES 

The  program  PBPEWAVE.M  displays  the  reflection  and  transmission  properties  of 
a  pipe  at  an  abrupt  change  in  diameter.  This  program  supports  instruction  in  concepts 
developed  in  Chapter  10,  Section  5  of  KFCS. 
A.         CONCEPTS 

The  general  equation  for  sound  power  reflection  coefficient  is  given  by 

(R0-pQc/S)2+X0i 

Rn= (11.1) 

(Vp0c/S)2+X02 

where  R^  and  Xq  are  resistance  and  reactance  of  the  termination. 

In  the  case  of  an  infinite  pipe  the  terminating  acoustic  impedance  is  purely 
resistive,  p0cS2,  simplifying  the  equation 

**iw  (u'2) 

The  power  transmission  coefficient  can  be  found  in  a  similar  way  and  is 

4SX 


T  - 


(VS2> 


(11.3) 


65 


The  power  transmitted  depends  only  on  the  cross  sectional  areas  of  the  two  parts 
of  the  pipe.  Therefore  the  power  transmission  for  any  pipe  in  any  medium  at  any  frequency 
can  be  found  if  the  ratio  of  areas  is  known. 


B. 


PROGRAM  DESCRIPTION 


The  program  PEPEWAVE.M  has  no  user  inputs  and  plots  the  reflection  and 
transmission  coefficients  for  an  infinite  pipe  with  a  range  of  cross  sectional  area  ratios. 
Figure  11.1  is  a  plot  of  the  reflection  and  transmission  coefficients.  Several  things  can  be 
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Figure  11.1  Reflection  and  Transmission  Coefficients  in  a  Pipe 
seen  from  this  figure.  As  the  two  sections  of  pipe  get  closer  to  the  same  size  the 
transmission  coefficient  becomes  one.  This  is  expected  as  it  is  one  long  pipe  of  constant 
cross  section  and  the  acoustic  impedance  remains  unchanged.  As  the  ratio  of  areas  tends 
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towards  zero,  the  reflection  coefficient  becomes  one,  which  is  the  same  as  a  pipe  being 

terminated  in  a  rigid  cap.  Notice  the  sum  of  reflection  and  transmission  coefficients  is 

always  one,  this  is  expected  as  no  loss  mechanisms  are  incorporated  in  the  model. 

C.         ALGORITHM  pipewave.m 

function  pipewave 

%  PIPEWAVE  computes  and  displays  reflection  and  transmission  coefficients  for  waves 

%  in  a  long  pipe  which  encounter  a  reduction  in  the  pipe  diameter.  The  program  requires 

%  no  user  input  and  plots  for  a  range  of  pipe  diameters. 

A=0:.01:l; 

R=((l-A).-2)./((l+A).-2); 

T=(4.*A)./((1+A)A2); 

figure(l) 

plot(A,R>:') 

hold  on 

plot(A,T>-') 

hold  off 

ylabel('Coefficient  Value') 

xlabel('Ratio  of  Second  Area  to  First  Area') 

legend('Reflection  Coeff,  'Transmission  Coeff) 
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Xn.  ACOUSTIC  FILTERS 

A.         LOW-PASS  FILTERS 
1.  Concepts 

A  short  section  of  large  diameter  pipe  of  length  L  and  area  S^ira2  inserted  in  a 
smaller  pipe  of  area  S  acts  as  a  low  pass  filter.  When  kL«l  and  ka«l,  the  reactance  is 

Substituting  into  the  equation  for  transmission  coefficient  of  a  branch 


71=- 


"b      "xi> 


Pnc  o  „  „  (I2-2) 

2S      b         b 


where  R^  and  Xb  are  the  input  resistance  and  reactance  of  the  branch  gives 

7,- ! 

l+(^)2  °23) 

2  S 

for  values  of  kL«l . 

A  more  complicated  analysis  results  in  an  equation  good  for  all  values  of  kL 
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Tn= 4- 

4cos2£Z,+(-+-)2sin2*Z,  <12-4) 

S    Sx 


Both  of  these  equations  are  used  in  the  program  LOPASSAC.M. 
2.  Program  Description 

LOPASSAC.M  requires  arguments  for  diameter  of  the  pipe  and  the  filter  section. 
Both  of  these  values  should  be  in  centimeters. 

The  program  computes  the  estimated  values  of  Tn  and  the  more  exact  solution  and 
plots  them  on  the  same  graph.  Figure  12. 1  is  a  plot  of  the  transmission  coefficient  for  a 
pipe  2  cm  in  diameter  and  a  filter  section  4  cm  in  diameter.  The  coefficients  are  plotted 
against  values  of  kL  from  0  to  1.  Notice  the  estimated  values  of  transmission  drop  off 
continuously  as  kL  increase.  However  the  more  exact  solution  shows  the  transmission 
coefficient  reaches  a  minimum  when  kL  is  tt/2,  or  one-quarter  wavelength.  It  then 
increases  to  a  maximum  of  1  when  kL  is  equal  to  a  half  wavelength.  This  pattern  continues 
for  values  of  kL  greater  than  1 .  One  other  interesting  aspect  of  this  program  is  that  it  can 
be  used  to  determine  the  transmission  coefficient  of  a  pipe  with  a  restriction.  The 
estimated  transmission  coefficient  will  only  be  good  for  very  small  values  of  kL  but  the 
coefficient  determined  by  the  more  exact  solution  is  correct  over  a  larger  range  of 
frequencies  as  long  as  ka«l . 
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Transmission  Coefficient  for  5  cm  Filter  in  3  cm  Pipe 
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Figure  12.1  Pipe  with  diameter  of  2cm  and  filter  section  4cm  in 
diameter. 

3.  Algorithm  lopassac.m 

function  lopassac(diam_pipe,diam_filt) 

%  LOP  AS  SAC  computes  and  plots  the  transmission  characteristics  of  a  low  pass  acoustic 

%  filter  in  a  pipe.  "diam_pipe"  and  "diam_filt"  can  be  any  units  the  user  wishes  as  long 

%  as  they  are  consistent. 

%  This  is  user  protection  portion  of  program. 

ifnargin~=2 

error('LOPASSAC  requires  two  input  arguments. ') 
end 
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if(min([diam_pipe,diam_filt])<=0) 

errorC  Arguments  for  LOP  ASS  AC  must  be  positive.') 
end 

%  This  begins  the  body  of  the  program. 
steps=200; 
domain=(l*pi); 
dom_step=domain/steps; 
plot_dom=domain/pi ; 
plot_dom_step=plot_dom/steps; 
kL=0:dom_step:domain; 
KL=0:plot_dom_step:plot_dom; 
Sl=((diam_filt)A2)*pi; 
S=((diam_pipe)A2)*pi; 

T_exac=4./((4.*(cos(kL))A2)+((((Sl/S)+(S/Sl))A2).*(sin(kL)).A2)); 
T_est=l./(l+((.5*(Sl/S).*kL).A2)); 
figure(l) 

plot(KL,T_exac,*-') 
hold  on 

plot(KL,T_est/:') 
hold  off 
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title( [Transmission  Coefficient  for  ',num2str(diam_fllt),... 

'  cm  Filter  in  ',num2str(diam_pipe),'  cm  Pipe']) 
xlabelCkL  Scaled  by  Pi') 
ylabel(Transmission  Coefficient') 
legend('Exact  Solution',  Long  Wavelength  Estimate') 

B.         HIGH-PASS  FILTERS 

1.         Concepts 

A  small  branch  of  length  L  and  radius  a,  in  the  wall  of  a  pipe,  will  act  as  a  high- 
pass  filter.  If  kL«l  and  ka«l,  the  branch  impedance  is 


p0ck2     PfX'w 
4rc         na2 


Substituting  this  back  into  equation  (12.2)  will  give  the  transmission  coefficient.  The 
ratio  of  branch  resistance  to  reactance  is 


R>    ka1 

(12.6) 


Xb    AU 


where  L'  is  the  effective  length  and  is  equal  to  L  +  1.5a.  If  ka2«l,  we  can  assume  the 
resistance  is  negligible  and  the  transmission  coefficient  is 
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r.— 1 


l+(J^L)2  (12-7) 


The  cutoff  frequency,  defined  as  the  frequency  where  T„  =  0.5,  is  given  by 


r     a2c 

f.-^j  (12.8) 


2.         Program  Description 

HCPASSAC.M  computes  the  estimated  transmission  coefficient  for  frequencies 
up  to  2500  Hz.  The  user  is  required  to  provide  four  input  arguments.  Branch  diameter 
and  length  and  pipe  diameter  should  be  in  centimeters  and  sound  speed  should  be  in 
meters/second.  The  program  plots  both  the  more  exact  transmission  coefficient  and  the 
estimated  coefficient. 

Figure  12.2  is  a  plot  for  a  6  cm  diameter  pipe  in  air  with  an  opening  3. 1  cm  in 
diameter  and  a  length  of  6  mm.  Looking  at  equation  (12.7)  attenuation  can  be  increased 
by  increasing  the  diameter  of  the  opening  or  by  decreasing  the  length  of  the  opening. 
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Highpass  Acoustic  Filter 
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Figure  12.2  High-pass  filter. 
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3.         Algorithm  hipassacm 

ftinctionhipassac(hole_diam,hole_length,pipe_diam,sndspd) 

%  HIPASSAC  computes  and  plots  the  transmission  characteristics  of  a  high  pass  acoustic 

%filter  in  a  pipe.  Hole  diameter,  hole  length  and  pipe  diameter  should  all  be  in  units  of 

%  centimeters.  Sound  speed  should  be  in  meters/second. 

%  This  is  user  protection  portion  of  program. 

ifnargin—4 

error('HIPASSAC  requires  four  input  arguments.') 
end 
if  (min([hole_diam,pipe_diam])<=0) 
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error('Hole  diameter  and  pipe  diameter  must  be  positive.') 
end 
if(hole_length<0) 

error('Hole  length  must  be  zero  or  greater.') 
end 
if(sndspd<=0) 

error  ('Sound  speed  must  be  positive.') 
end 

%  This  begins  the  body  of  the  program. 
freqs=logspace(  1 ,4, 1 00); 
w=2*pi.*freqs; 
k=2*pi.*freqs./sndspd; 
S=((pipe_diam/(2*  100))A2)*pi; 
hole_length=hole_length/100; 
a=hole_diam/(2*100); 
L_effHiole_length+(1.5*a); 
rho=1.2; 

R=rho*sndspd.*(kA2)./(4*pi); 
X=rho*L_eff.*w./(pi*(aA2)); 
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T_est=l./(l+(((pi*(aA2))./(2*S*L_eff.*k)).A2)); 

T^(R.A2)+(X  A2))./((((rho*sndspd/(2*S))  +  R).A2)+(X.A2)); 

figure(l) 

semilogx(freqs,T_est,'-') 

hold  on 

semilogx(freqs,T,':') 

hold  off 

xlabel(Trequency  (Hz)') 

ylabel(Transmission  Coefficient') 

title(TIighpass  Acoustic  Filter') 

legendCLow  Freq  Limit',  *Exact  Coefficient') 

C.         BAND-STOP  FILTERS 

1.         Concepts 

A  Helmholtz  resonator  attached  to  a  pipe  acts  as  a  band-stop  filter.  The  resonator 
does  not  absorb  sound  so  resistance  is  zero  and  the  reactance  is 


v       ,  ixiL      c    >. 


where  V  is  the  volume  of  the  resonator  and  Sr  is  the  area  of  the  resonator  opening. 
Substituting  into  equation  12.2  the  transmission  coefficient  is 
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T  = 


U c- 


wZ/     c2  (12-10) 

4S2(---)2 
Sr     uV 


The  center  frequency  of  the  band-stop  filter  is  given  by  the  resonant  frequency  of  the 
resonator 


f,-c 


S 

(12.11) 


2n\  L'V 

2.         Program  Description 

The  program  BDSTOPAC.M  requires  four  input  arguments.  Neck  diameter  and 
length,  and  pipe  diameter  should  be  in  centimeters  and  volume  should  be  cubic 
centimeters.  Figure  12.3  is  a  plot  of  a  band-stop  filter  with  an  opening  3.1  cm  in 
diameter,  a  neck  length  of  6  mm,  a  volume  of  1 120  cc,  and  a  pipe  diameter  of  6  cm. 
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Band-stop  Acoustic  Filter 
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Figure  12.3  Band-stop  filter. 


3.  Algorithm  bdstopac.m 

function  bdpassac(hole_diam,hole_length,res_vol,pipe_cham) 

%  BDSTOPAC.M  computes  and  plots  the  transmission  characteristics  of  a  band  pass 

%  acoustic  filter  in  a  pipe.  Hole  diameter,  hole  length  and  pipe  diameter  should  all  be  in 

%  units  of  centimeters.  Resonator  volume  should  be  in  units  of  cubic  centimeters. 

%  This  is  user  protection  portion  of  program. 

if  nargin~=4 


errorCBDSTOPAC  requires  three  input  arguments') 


end 


if(min([hole_diam,pipe_diam,res_vol])<=<)) 
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errorCHole  diameter,  pipe  diameter  and  resonator  volume  must  be  positive.') 
end 
if(hole_length<0) 

errorCHole  length  must  be  zero  or  greater.') 
end 

%  This  begins  the  body  of  the  program. 
freqs=l:2500; 
omega=2*pi.  *freqs; 
sndspd=343; 

S=((pipe_diam/(2*  100))A2)*pi; 
hole_length=hole_length/l  00; 
a=hole_diam/(2*100); 
L_effHiole_length+(  1 .  7*a); 
Sb=pi*(aA2); 
vol=res_vol/(  1  e6); 
T=l./(l+((sndspdA2)./(4*(SA2).*((((omega.*L_eff)./Sb)-... 

((sndspdA2)./(omega.*vol)))A2)))); 
figure(l) 

semilogxCfreqs,!,'^) 
xlabel(Trequency  (Hz)') 


80 


ylabel(Transmission  Coefficient') 
titleCBand-stop  Acoustic  Filter') 
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Xffl.  RECEIVER  OPERATOR  CHARACTERISTICS 

A.  CONCEPTS 

The  program  ROC.M  and  its  several  supporting  programs  are  designed  to  allow 
the  user  to  determine  Receiver  Operator  Characteristics  curves.  This  program  supports 
instruction  in  concepts  discussed  in  Chapter  1 1  of  KFCS. 

B.  PROGRAM  DESCRIPTION 

This  program  actually  consists  of  two  main  programs  and  5  smaller  support 
programs.  Figure  1  is  the  user  interface  that  is  generated  by  the  program  ROC.M. 


RECEIVER  OPERATING  CHARACTERISTICS 

USER  CONTROLS 

jj  J  jj 

200  Hz  1100  2  kHz 

Frequency 


LEARNING 

While  Noise 


1  kHz  Tonal 
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New  Lab 
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Tonal  Amp  (S/N) 


Next  Signal 


Repeat 


Signal  Present? 


YES 


NO 


Figure  13.1  ROC  User  Interface. 


When  the  program  is  started  the  pushbuttons  "YES",  "NO",  "NEXT  SIGNAL' 
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and  "REPEAT'  are  not  yet  needed  and  cannot  be  seen  by  the  user.  Throughout  the  run 
when  a  pushbutton  is  not  needed  it  is  not  visible  to  the  user. 

The  slide  bars  "Frequency"  and  "Tonal  Amp  (S/N)"  allow  the  user  to  adjust  the 
frequency  and  amplitude  of  the  tonal.  The  value  can  be  changed  by  moving  the  slide  bar 
or  by  depressing  the  button  at  either  end  of  the  slide  bar.  The  frequency  and  amplitude 
are  displayed  below  the  slide  bar.  These  sliders  can  be  moved  at  anytime  however,  once 
a  run  is  started  the  frequency  and  amplitude  of  the  tonal  cannot  change  until  the  run  is 
complete. 

For  training  purposes,  "White  Noise"  and  the  "1  kHz  Tonal"  pushbuttons  allow 
the  user  to  listen  to  the  noise  or  tonal  separately.  The  "Tonal  +  Noise"  pushbutton  allows 
the  user  to  listen  to  the  two  signals  combined  at  equal  amplitudes.  These  three  buttons 
can  be  pressed  at  anytime  during  a  run  but  the  frequency  of  the  tone  will  always  be  1 
kHz. 

"New  Lab"  starts  a  run  and  makes  the  "Next  Signal"  burton  visible.  Currently  the 
program  creates  ten  signals  per  run,  to  make  the  data  more  statistically  significant  the 
program  variable  "trials"  in  NEWLAB.M  can  be  changed  to  any  number  the  user  wishes. 
All  signals  contain  random  noise  and  the  program  makes  a  random  yes/no  decision  for 
each  signal  to  determine  whether  it  will  contain  a  tonal.  If  the  decision  is  yes  the  program 
adds  a  tonal  at  the  frequency  and  amplitude  specified  by  the  user  to  the  signal. 

When  the  "Next  Signal"  is  depressed  the  first  signal  will  be  played  for  four 
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seconds.  The  "YES",  "NO",  and  "Repeat"  buttons  then  show  up.  The  user  may  elect  to 
hear  it  again  by  depressing  "Repeat".  After  the  user  depresses  the  "YES"  or  "NO"  button 
the  next  signal  can  be  heard  by  again  depressing  the  "Next  Signal"  button.  After  the  run 
is  complete  the  buttons  will  again  become  invisible  and  the  "New  Lab"  will  be  visible. 

After  completion  the  results  of  the  run  will  be  displayed.  The  number  of  signals, 
number  with  a  tone  present,  number  correctly  identified,  number  with  no  tone  that  were 
incorrectly  identified,  a  probability  of  detection,  and  a  probability  of  false  alarm  will  all 
be  displayed. 

ROC  curves  can  now  be  generated  using  information  from  this  lab.  However,  it  is 
important  to  realize  several  things.  First,  it  takes  several  points  to  accurately  define  any 
single  curve.  To  do  this  it  is  best  to  keep  the  S/N  ratio  constant  during  several  trials  but 
change  the  criterion  for  answering  "YES".  The  first  time  the  criterion  might  be  to  answer 
yes  if  there  is  even  the  slightest  possibility  there  is  a  signal.  The  next  criterion  might  be 
to  say  yes  only  if  you're  positive  there  is  a  signal.  This  will  give  several  points  along  the 
same  curve.  Other  curves  can  be  determined  by  changing  the  S/N  ratio  and  using  the 
same  yes/no  criteria. 
C.        ALGORITHMS  FOR  ROC  LAB. 

1.  Algorithm  roc.m 

%function  roclab 
%  This  program  is  the  first  to  be  called  for  the  ROC  lab.  It  generates  the  user  interface  . 
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%  It  requires  no  user  protection  because  the  interface  itself  only  allows  certain  inputs  to 

%  be  made. 

clear 

global  trial 

global  signals 

global  trials 

global  resp 

noisel=randn(l,((2A13)*4)-4); 

freql=500; 

arg_step=(2*pi*rreql*4)/(((2A13)*4)-l); 

sin_arg=0:arg_step:((2*pi*freql*4)-arg_step); 

sigl=sin(sin_arg); 

sig2=sigl(l,l:(length(sigl)-3)); 

mix_sig=noise  1 +sig2; 

rocwin^gureCunitsVnormalVpositiori^^  .25  .5  .5],... 

,name',,Receiver  Operator  Characteristics  Lab',... 

'colorVwhiteVnumbertitleVofF); 
axis  off 
textCunits'/normalizedVposition',  [  .5.95],... 

TiorizontalalignmentVcenterVstring',... 
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■RECEIVER  OPERATING  CHARACTERISTICS7color',*black'); 

textOumtsVnormalizedVpositioriX  1  .85],... 

'horizontalalignmentVcenterVstring',. . . 

•LEARNING'/color'^lack'); 
text('unitsVnormalizedVposition',[.7  .85],... 

'horizontalalignmentVcenterVstring',. . . 

USER  CONTROLSVcolor'^lack'); 
only_noise^control('style7pushbuttori,'units',. . . 

'normalized','positiori,[.  1 1  .65  .225  .1],... 

rbackgroundcolorVmagenta',,foregroundcolor','black',.  ■ . 

lhorizontalalignment','centerVstring,,'White  Noise',. . . 

'callback','sound(noise  1 );'); 
onlytone^^icontrolCstyle'/pushbuttonVunits', . .  . 

'normalizedVposition',[.  1 1  .5  .225  .1],... 

'backgroundcolor','magenta,,,foregroundcolor',fblack', . . . 

'horizontalalignmentVcenter','stringVl  kHz  Tonal',.-. 

'callback','sound(sig2);'); 
tone_and_noise=uicontrol('style','pushbuttonVunits',. . . 

'noraialized','position',[.  1 1  .35  .225  .1],... 
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,backgroundcolo^Vmagenta',,fo^egroundcolor',,black,,. . . 

,horizontalalignment','center,,,string,,Tonal  +  Noise',... 

'callbackVsound(mix_sig);'); 
new_lab=uicontrol('styleVpushbuttonVunits',. . . 

'noraialized','position',[.  1 1  0.2  .225  .1],... 

,backgroundcolorVmagentaVforegroundcolor',rblack',. . . 

'horizontalalignmentVcenterVstring',. . . 

TSfew  Lab'/callback'/newlab'); 
user_freq^icontrol('styleVslider7units7normalized',. . . 

'position',[45  .7  .45  .05],rbackgroundcolor,,'magenta',... 

Toregroundcolor\'black7min,,200,,max',2000,... 

Value'JlOO/callbackVfreqdisp'); 
u_f=num2str(get(user_fTeq,Value')); 
text('imits7normalized7position',[.7.675],... 

Tiorizontalalignment','centerVstring,,u_f,,color,,,black'); 
text('units7normalized7positiori,[.5  .675],... 

Tiorizontalalignment'/rightVstring',. . . 

'200  Hz7color*,'black'); 
text('imits7normalized7position',[.9  .675],... 

TiorizontalalignmentVleftVstring1,. . . 
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'2  kHzVcolor^lack'); 
text('unitsVnormalizedVposition',[.7  .6],... 

'horizontalalignmentVcenterVstring',. . . 

TrequencyVcolor','black,); 
tonal_amp=uicontrol(,styleVslider,,'unitsVnormalized',. . . 

'position', [.4 5  .475  .45  .05],'backgroiindcolor',,magenta',. 

'foregroundcolor','blackVmin',0,'max',  1 ,. . . 

'value',0.5,'callback','ampdisp'); 
t_a=num2str(get(tonal_amp,'value')); 
text('units ','nonnalizedVposition',  [.  7.4],... 

'horizontalalignmentVcenter,,'string',t_a,'color','black'); 
text('units','normalizedVposition',[.45  .4],.. . 

'horizontalalignment','right','string',. . . 

'OVcolor'^lack'); 
text('unitsVnormalized','position',[.95  .4],... 

'horizontalalignment','left ','string',. . . 

'lVcolor'^lack'); 
textCunits'/normalized'/position'^^  .325],... 

TiorizontalalignmentVcenterVstring',. . . 

Tonal  Amp  (S/Ny/color'^lack'); 


89 


text('unitsVnormalized','position',[0. 15  0],... 

'horizontalalignmentVcenterVstring',. . . 

'Signal  Prasei^'/color'/black*); 
s  ig_yes=uicontrol(,styleVpushbuttonVunits',. . . 

'normalizedVposition*^^  0.05  .2  .1],... 

TbackgroundcolorVmagentaVfo^egroundcolor','black,,. 

Tiorizontalalignment'/centerVstring',. . . 

,YES,;callback,,,resp=l;results;,,Visible,,,ofr); 

sig_no=uicontrol(,styIeVpushbuttonVunits,v . . 

'normalizedVposition',[.7  0.05  .2  .1],... 

rbackgroundcolorVmagentaVforegroundcolor','black',. 

'horizontalalignment'/centerVstring',. . . 

TSfO,;callback,;resp=2;results;,,Visible',,off); 
next_sig=uicontrol('styleVpushbuttonVunits',. . . 

'normalized','position',[.45  0.2  .2  .1],... 

fbackgroundcolorVmagentaVforegroundcolor','blackV 

*horizontalalignment',lcenter','string',. . . 

Tvfext  SignalVcallbackVnextsig;',..- 

VisibleVofF); 


90 


repeat_sig=uicontrol(,style,,'pushbutton',lunits',. . . 

'normalizedVposition',[.7  0.2  .2  .1],... 

rbackgroundcolor7magenta7foregroundcolor,,rblack',. . . 

'horizontalalignmentVcenterVstring',. . . 

"RepeatVcallbackVsoundC  signals(trial, :));',... 

Visible'/off); 

2.         Algorithm  newlab.m 
%function  newlab 

%  This  program  generates  the  set  of  signals  to  be  used  during  a  run  for  the  roc  lab. 
global  trials 
global  signals 
global  trial 
global  tonalcoeff 

%  The  following  section  sets  defaults  for  variables  used  in  program 
trial=l: 


frequency=get(user_freq,'value'); 


in  Hz 


%  frequency  of  tonal 


noise  duration=4: 


trials=10; 


sig_mag=get(tonal_amp, 'value'); 


%  duration  of  sound  in  seconds 


%  relative  magnitude  of  tonal 
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noise_length=(2  A 1 3  )*noise_duration;  %  this  is  total  number  of  samples 

%  that  will  be  used,  2A13  per  second 
%  determined  by  sampling  rate  of 
%  sound  function  which  is  8192  Hz 

noise=Tandn(trials,noise_length);  %  Generates  white  noise  signal 

%  The  following  section  creates  the  tonal  and  determines  which  signals  will  have 

%  the  tonal  present. 

a=(2*pi*user _freq*noise_duration)/(size(noise,2)-l);  %  finds  the  step  size  for  tonal 

%  generation 

freq_count=(2*pi*user_freq*noise_duration); 

sig=0:a:freq_count;  %  generates  all  pts  for  argument  of  sin 

tonal=2*sig_mag*sin(sig); 

rand('seed,,sum(  1 00*clock)); 

tonal_coeff=round(rand(trials,  1 )); 

tone_matrix=(tonal_coeff*tonal); 

signals=noise+tone_matrix; 

setCnextsig/visible'/on*); 

set(new_lab,Visible','ofF); 

refresh(rocwin); 

3.         Algorithm  nextsig.m 
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%function  nextsig 

%  This  program  produces  the  tone  to  begin  each  signal  for  the  roc  lab. 

sound(signals(trial,:)); 

set(next_sig,'visibleVoff); 

set(repeat_sig,'visible7on'); 

set(sig__yes,'visibleVon'); 

set(sig_no,'visibleVon'); 

refresh(rocwin); 

4.  Algorithm  freqdisp.m 

%  This  is  a  script  file  used  to  display  the  current  frequency. 

prev_f=u_f; 

u_f=num2str(get(user_rreq,'value')); 

text('units7nomialized7positiori,[.7  .675],... 

'horizontalalignment7center7srring',prev_f,. . . 

'color'/white'); 
text(,units',,normalized7position,,[.7  .675],... 

Tiorizontalalignment7center','string,,u_f,lcolor',fblack'); 

5.  Algorithm  ampdisp.m 

%  This  is  a  script  file  used  by  thr  roc  lab  to  display  the  S/N  ratio. 
prev_ta=t_a; 
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t_a=num2str(get(tonal_amp,'value')); 
text('units7normalized7position',[.7  .4],... 

'horizontalalignmentVcenterVstring^prevta,. . . 

'color'/white'); 
text('iinits7normalized7positiori,  [.  7  .4],. . . 

'horizontalalignmentVcenterVstring,,t_a,'color',,black'); 

6.         Algorithm  results.m 
%function  results 

%  This  function  is  used  by  the  roc  lab  to  compute  results  at  the  end  of  a  run  of 
%  frequencies, 
global  resp 
global  tonalcoeff 
global  trial 
global  trials 
if  trial=l 

detect=0; 

sig_pres=0; 

false_alarm=0; 

sig_notpres=0; 
end 
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if  resp==  1  %  Answer  is  yes 

sig_there=l; 
elseif  resp=2  %  Answer  is  no 

sig_there=0; 
end 

if  tonal_coefT(trial)=  1 
ifsig_there=l 

detect=detect+l; 
sig_pres=sig_pres+ 1 ; 
elseif  sig_there=0 

sig_pres=sig_pres+ 1 ; 
end 
elseif  tonal_coeff(trial)=0 
if  sig_there=0 

sig_notpres=sig_notpres+ 1 ; 
elseif  sig_there=  1 

sig_notpres=sig_notpres+ 1 ; 
false_alarm=false_alarm+ 1 ; 
end 
end 
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iftrial=trials 

fprintf('%3. Of  trials  were  run.\n',trials) 

fprintf('%3. Of  trials  had  a  signal  present.  \n',sig_pres) 

fprintf('You  correctly  identified  %3.0f  of  these  signals.\n',detect) 

fprintf('In  addition  %3.0f  of  the  trials  were  incorrectly\n',false_alarm) 

fprintfC  identified  as  having  a  signal  present.W) 

fprintf(This  results  in  a  Pd  of  %1.2f  \n',(detect/sig_pres)) 

fprintf('And  a  Pfa  of  %1.2f  \n',(false_alarm/sig_notpres)) 

se^newlab/visible'/on'); 

set(sig_yes,'visibleVoff); 

set(sig_no,'visible','ofF); 

set(repeat_sig,Visible7off); 

else 

setCnextsig/visible'/on'); 

set(repeat_sig,'visible','off); 

trial=trial+l; 

end 

refresh(rocwin); 
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APPENDIX.  MODIFIED  POLAR  PLOT  PROGRAM  LISTING 

function  hpol  =  polar(theta,rho,line_style) 

%POLAR        Polar  coordinate  plot. 

%         POLAR(THETA,  RHO)  makes  a  plot  using  polar  coordinates  of 

%         the  angle  THETA,  in  radians,  versus  the  radius  RHO. 

%         POLAR(THETA,RHO,S)  uses  the  linestyle  specified  in  string  S. 

%         See  PLOT  for  a  description  of  legal  linestyles. 

% 

%         See  also  PLOT,  LOGLOG,  SEMILOGX,  SEMILOG Y. 

%         Copyright  (c)  1984-94  by  The  MathWorks,  Inc. 

%  This  program  has  been  modified  to  be  used  with  the  program  BKNARRAY.M 

%  by  Kevin  Melody. 

global  radius 

global  r3db 

if  nargin<  1 

errorCRequires  2  or  3  input  arguments.') 
elseif  nargin  =  2 
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if  isstr(rho) 

linestyle  =  rho; 
rho  =  theta; 
[mr,nr]  =  size(rho); 
ifmr=  1 

theta  =  1  :nr; 
else 

th  =  (l:mr)'; 

theta  =  th(:,ones(l,nr)); 
end 


else 


linestyle  =  'auto'; 


end 
elseifnargin  =  1 

linestyle  =  'auto'; 
rho  =  theta; 
[mr,nr]  =  size(rho); 
if  mr=  1 

theta  =  1  :nr; 
else 
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th  =  (l:mr)'; 

theta  =  th(:,ones(l,nr)); 

end 
end 
if  isstr(theta)  |  isstr(rho) 

error('Input  arguments  must  be  numeric.'); 
end 
if  any(size(theta)  ~=  size(rho)) 

error('THETA  and  RHO  must  be  the  same  size.'); 
end 

%  get  hold  state 

cax  =  newplot; 

next  =  lower(get(cax,'NextPlot')); 

holdstate  =  ishold; 

%  get  x-axis  text  color  so  grid  is  in  same  color 
tc  =  get(cax,'xcolor'); 

%  Hold  on  to  current  Text  defaults,  reset  them  to  the 
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%  Axes'  font  attributes  so  tick  marks  use  them. 

fAngle  =  get(cax,  DefaultTextFontAngle'); 

fName   =  get(cax,  DefaultTextFontName'); 

fSize   =  get(cax,  DefaultTextFontSize'); 

fWeight  =  get(cax,  DefaultTextFontWeighf); 

set(cax,  DefaultTextFontAngle',  get(cax,  TontAngle'), ... 

'DefaultTextFontName',    get(cax,  TontName'), ... 

DefaultTextFontSize',    get(cax,  TontSize'), ... 

DefaultTextFontWeight',  get(cax,  TontWeight') ) 

%  only  do  grids  if  hold  is  off 
if~hold_state 

%  make  a  radial  grid 
hold  on; 

hhh=plot([0  max(theta(:))],[0  max(abs(rho( :)))]); 
v  =  [get(cax,'xlim')  get(cax,'ylim')]; 
ticks  =  length(get(cax,'ytick')); 
delete(hhh); 
rmin=0: 
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rmax=radius; 

%  define  a  circle 

th  =  0:pi/50:2*pi; 
xunit  =  cos(th); 
yunit  =  sin(th); 
%  now  really  force  points  on  x/y  axes  to  lie  on  them  exactly 
inds  =  [l:(length(th)-l)/4:length(th)]; 
xunits(inds(2:2:4))  =  zeros(2,l); 
yunits(inds(l:2:5))  =  zeros(3,l); 

%rinc  =  (rmax-rmin)/rticks; 
for  i=rmax 

plot(xunit*i,yunit*i,'-Vcolor,,tc,,liriewidth',  1 ); 

%text(l,i+rinc/20,['  '  num2str(i)],'verticalalignment','bottom, ); 
end 
for  i=r3db 

plot(xunit*i,yunit*i,':  Vcolor',tc,'linewidth',  1 ); 

%text(l,i+rinc/20,['  '  num2str(i)],Verticalalignment',,bottom, ); 
end 
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%  plot  spokes 

th  =  (l:2)*2*pi/4; 

est  =  cos(th);  snt  =  sin(th); 

cs  =  [-est;  est]; 

sn  =  [-snt;  snt]; 

plot(nnax*cs,rmax*sn,'-7color',tc,'linewidth',  1 ); 

%  annotate  spokes  in  degrees 
rt=  1.2*rmax; 
for  i  =  1  :max(size(th)) 

text(rt*cst(i),rt*snt(i),int2str(i*90),,horizontalalignmentVcenter' ); 
if  i  =  max(size(th)) 

loc  =  int2str(0); 
else 

loc  =  int2str(180+i*90); 
end 

text(-rt*cst(i),-rt*snt(i),loc,'horizontalalignmentVcenter'); 
end 

%  set  viewto  2-D 
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view(0,90); 
%  set  axis  limits 

axis(rmax*[-l  1-1.1  1.1]); 
end 

%  Reset  defaults. 

set(cax,  'DefaultTextFontAngle',  fAngle  , ... 

T)efaultTextFontName',    fName , ... 

TtefaultTextFontSize',    fSize, ... 

TDefaultTextFontWeighf,  fWeight ); 

%  transform  data  to  Cartesian  coordinates. 
xx  =  rho.*cos(theta); 
yy  =  rho.*sin(theta); 

%  plot  data  on  top  of  grid 
if  strcmrXlinestyle/auto') 

q  -  plot(xx,yy); 
else 

q  =  plot(xx,yy,line_style); 


103 


end 

if  nargout  >  0 

hpol  =  q; 
end 
if~hold_state 

axis('equar);axis('off); 
end 

%  reset  hold  state 

if -holdstate,  set(cax,TSIextPlot',next);  end 
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